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ABSTRACT 

A micron n + — n — n + silicon diode is simulated via the hydrodynamic model for carrier 
transport. The numerical algorithms employed are for the non-steady case, and a limiting 
process is used to reach steady state. The novelty of our simulation lies in the shock capturing 
algorithms employed, and indeed shocks, or very rapid transition regimes, are observed 
in the transient case for the coupled system, consisting of the potential equation and the 
conservation equations describing charge, momentum, and energy transfer for the electron 
carriers. These algorithms, termed essentially non-oscillatory, have been successfully applied 
in other contests to models the flow in gas dynamics, magnetohydrodynamics and other 
physical situations involving the conservation laws of fluid mechanics. The method here is 
first order in time, but the use of small time steps allows for good accuracy. Runge-Kutta 
methods allow one to achieve higher accuracy in time if desired. The spatial accuracy is of 
high order in regions of smoothness. 
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1. INTRODUCTION 


A micron n + — n — n + silicon diode is simulated via the hydrodynamic model for 
carrier transport. 'The numerical algorithms employed are for the non-steady case, 
and a limiting process is used to reach steady state. The novelty of our simulation 
lies in the shock capturing algorithms employed, and indeed shocks, or very rapid 
transition regimes, are observed in the transient case for the coupled system, con- 
sisting of the potential equation and the conservation equations describing charge, 
momentum, and energy transfer for the electron carriers. These algorithms, termed 
essentially non-oscillatory, have been successfully applied in other contexts to model 
the flow in gas dynamics, magnetohydrodynamics and other physical situations in- 
volving the conservation laws of fluid mechanics. The method here is first order in 
time, but the use of small time steps allows for good accuracy. Runge-Kutta meth- 
ods allow one to achieve higher accuracy in time if desired. The spatial accuracy is 
of high order in regions of smoothness. 
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The conservation laws for charge, momentum, and energy for each carrier can 


be obtained by taking the first three moments of the Boltzmann equation [see 
Blotekjaer [3], Rudan and Odeh [23], and Cercignani [4]]. This will result in a 
system of five equations with fourteen variables in three dimensional space. These 
variables are defined by particle density n, velocity v, internal energy e/, heat flux 
q, and symmetric pressure tensor p. In order to close the system, we assume that p 
and q are defined via (cf. Gardner, Jerome and Rose [9]), 


hnT _ -k „ 

Pij — ®ijr 

m m 


where kb is the Boltzmann constant, m is the mass of one carrier, and T is the 
carrier temperature assumed to be scalar. The moment system also includes aver- 
age collision terms, denoted C n ,C p , and C w , respectively, which describe the rate 
of change of mass, momentum and energy. They will account for recombination 
regeneration processes, electron-electron and electron-lattice collisions, and transfer 
of energy between electrons and lattice. Their explicit form, for the application 
considered here, is given by 

C n = 0, 


Cp = 


-P 


C w = 


~(w - jnkbTo) 


(P := carrier momentum density) 

(w := carrier energy density) 


Here we use the following empirical relations [see [9]] for relaxation times and the 


thermoconductivity constant: 


__ mpoTo _ _ SpokbTTo “ 1 r p ^ _ 3/x 0 fc|7o 

Tp ~ eT ’ Tw ~ 2ev 2 s {T + T 0 ) + 2 ’ ^ “ 2e 
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where To is the temperature of the lattice (possibly a function of space), hq is the 
electron mobility, e the charge of an electron, and v 3 is the saturation velocity. We 
shall discuss the foundations of these relations below. We note here, for consistency, 
that e/ is proportional to T, and in the following section we shall define P and it;, 
and describe the equations of the hydrodynamic system and the system boundary 
conditions. 

The recent literature, pertaining to the device model simulated here, includes, 
in addition to [23] and [9], the studies by Baccarani and Wardeman [1], by Odeh, 
Rudan, and White [19], and by Fischetti and Rudan [8]. In each case, velocity 
overshoot effects are simulated, which cannot be recovered from the standard drift- 
diffusion 'model. As is now widely understood, the effect is strongly dependent on 
the semiconductor, and is much more pronounced in gallium arsenide than silicon, 
for example. Since velocity overshoot is related to variations in termperature and 
electric field, the energy equation in the hydrodynamic model thus allows for the 
heat flux term q governed by the Wiedemann-Franz law [2], (note that kti is the 
conductivity in this law) and the coupling to the potential equation is essential in 
this case. Moreover, the momentum equation plays the role of a refined constitutive 
current relation. Rather than seeking modifications of the classical current relation 
by the incorporation of relevant physical effects beyond critical values of parameters 
such as i (see [1] for such values and Thomber [30] for a sophisticated discus- 
sion of such augmented relations), we find here that current is induced by variable 
concentration and momentum within the system. 
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The validity of the assumptions concerning the collision terms C p and C w , 
termed phenomenological approach in [7], has been discussed by Nougier et al. [18], 
where agreement within 10%, for n type silicon and p type germanium, has been 
demonstrated, relative to a direct simulation of the Boltzmann equation, for a wide 
range of electric fields and lattice temperatures. Further assumptions employed 
in the model, also used in [1], [23], [19], and [9], are the existence of a scalar 
quantity for the temperature tensor, constructed from the random velocity, and 
the parabolic representation for the energy, in terms of its kinetic and thermal 
components. Following Hess and Sah [13], and [1], we have assumed a constant 
diffusivity coefficient, which is consistent with the assumptions of constant effective 
mass, and scalar temperature and diffusivity. This has the effect, via the Einstein 
relation, of describing the electron mobility as an inverse temperature dependent 
function. 

Thus, in this approach, mobility does not directly depend upon velocity. This 
mobility relation explicitly dictates the expression for the temperature dependent 
momentum relaxation times expressed above, as well as the choice, exponent = 
— 1, in the Wiedemann-Franz law for the thermal conductivity. The temperature 
dependent energy relaxation time is deduced by using the empirical expression for 
the mobility, presented by Caughey and Thomas [5], so as to express drift velocity 
in terms of temperature. The resulting expression is therefore empirical. The 
steady-state character of both the resulting relaxation time expressions has been 
emphasized in [1]. Since we simultaneously obtain information for the transient 
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analysis, the latter results must be interpreted provisionally with respect to t p and 
T w . 

The earliest theoretical predictions of velocity overshoot appear to be due to 
Rees [20], for germanium, and to Ruch [22], for silicon and gallium arsenide. Fur- 
ther studies by Shur [26] and Maloney and Frey [16] gave a detailed qualitative 
understanding of the phenomenon; in particular, it was understood that the physi- 
cal principle underlying velocity overshoot is the finite time taken before the energy 
of the carriers reaches equilibrium with the new field, as they traverse a region of 
rapidly varying electric field. These ideas, and the relation to ballistic transport, 
have been discussed and critiqued by several authors (see Shur and Eastman [27] 
and Hess [12]). 

Increasing sophistication of the model followed. While [26] and [16] were based 
upon the quasi-diode, by 1982 Cook and Frey [6] had incorporated many essential 
features of the hydrodynamic model, coupled to the potential equation, in their 
simulations. However, they did make the assumptions not made here, of vanishing 
heat flux and negligible convective energy term. The former assumption, reducing 
order, substantially alters the mathematical character of the energy equation in ad- 
dition to its physical implications; these mathematical implications are manifested 
principally in two or more dimensions, however. Other attempts to incorporate 
energy and/or transport effects followed, based upon various models. Thus, the 
foundation for. the work of Tang [29] was the spectral method introduced by Strat- 
ten [28], where a spherical harmonic expansion of the distribution function leads 
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to momentum and energy relations similar to those subsequently derived for the 
hydrodynamic model.. Also, McAndrew, Singhal, and Heasell [17] employ a model 
involving the momentum equation, but not the energy equation. 

The models employed in [1], [23], 19], and [9] are virtually the same, with some 
variations due to mobility-temperature representations, and the joint connection 
with relaxation times and the Wiedemann- Franz law. The model in this paper is 
the same as that of [9]. The numerical methods in this paper, however, are quite 
different from those in previously reported work. Both the parametric regimes and 
the iterative differencing schemes used in previous work were chiefly keyed to the el- 
liptic or elliptic/parabolic classification of the underlying equations, in other words, 
to what might loosely be characterized as the subsonic transport regime, where con- 
centrations and velocities would not be expected to experience shock discontinuities. 
In steady-state, the relevant “soundspeed” is y/ngT/m, obtained by linearization of 
the density-momentum subsystem. If the density-momentum-energy subsystem is 
linearized, the traditional “soundspeed” of is obtained, but this involves 

the parabolic mode associated with the energy equation. The analysis of [9], in fact, 
shows that a damped Newton method can be justified by energy type estimates in 
the subsonic case. The shock capturing methods used herein have both the ad- 
vantage of simulating broad parameter ranges, including transonic and supersonic 
flow, as well as being readily extendible to two dimensions. Since, as noted by [1], 
velocity overshoot is maximally exploited when cold carriers are injected at high 
field, the simulation should have the capability of covering a broad temperature 
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range, which suggests regions of transonic flow and complex behavior. 

2. SYSTEM AND BOUNDARY CONDITIONS 

The moment equations are specified, via the summation convention in three 
dimensions, as 

d t (n) + d Xi (nvi) = C n 

dt(nvj) + d Xi (nv{Vj + Pij ) = nFj + Cj j th component 

dt(n(~v 2 + e/)) + d X{ (nvi(^v 2 + e/) + PijVj + q <) = nFiV { + C w 

where all of the terms have been defined except the forcing function F due to the 
electric field. 

The Poisson equation for the electrostatic potential, the first of the Maxwell 
equations, is given in a one carrier, electron system by 

eV 2 $ = e(n — no) 


where rip is the density of donors. Thus the forcing function F becomes: 

e — e 

F = — E = — V* , 
m m 

We define P = mnv and w = mnej + m = | nki,T 4- \mnv 2 . Let the vector 
u = (n, P, w) denote, respectively, density, momentum, and energy of electrons. We 
get the final form of our equations for electrons and electric potential as: 

d t (n) +d Xi (nvi) = C n 
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dt(Pj) + d Xi (PiVj + hnT ) = -enEj + C p 
d t (w ) + d Xi (vi(w + kf,nT )) = —env{Ei + C w + d Xi (Knd Xi T) 

eV 2 $ = e(n — no) 

We can write the equations in the compact form in ( x,y,z ) space as: 

u t + /i(«)x + f 2 (u) y + fz{u) x = C(u) + G(u, $) + (0,0, V(/cnVT)) 

eV 2 $ = e(n — no) 


In this paper then we concentrate on solving the following system: 


(la) 

(lb) 


e$ xx =e(n - n D ), E = -$ x , 
n t + (nv) x = 0, 


P t + ( Pv -f k h nT) x = -enE , 

r p 

. , w - lnk b T 0 

w t + (vw + vkf,nT) x = —tnvE 




with the following boundary conditions 


$( 0 ) = bLi n ^!M) 

e n, 

$(1) = ^^/n(^^-) + vbias(t) 
e n,j 

u z (0) = Ux(l-) = 0. 


This is an incompletely parabolic system, very similar to the Navier- Stokes 
equations of compressible flow. Boundary conditions of this type for that system 
have been discussed, e.g. in [10]. 
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3. NUMERICAL SCHEME 


The numerical scheme is based on the essentially non-oscillatory shock captur- 
ing algorithms developed in [11]. We use a sixth order stencil in sp&ce and forward 

U ^ *** U ^ 

Euler in time. Let u" = u(jAx,nAt) = u(xj,t„). Also let D_u” = , 

D + U l = Uj± KT l ^ D o «? = 


Then using the above conventions we can write our scheme as: 


(2a) <•£>+£_$’• = e(n" - n Dj ) 

(2b) E” = - D 0 

(2c) «" +I - uj - - /,•_*) + AiC; + A tGJ + AtD.(nD + T) 

We calculate the potential from th© density and then we calculate the electric 
field. We calculate the solution in the next time step from the electric field and the 
’’flux” calculated in the previous step. 

The “flux” term, /”_ t , represents the numerical flux flowing into the jth 

j ^ 

computational cell [see Lax [15]]. 

The finite difference scheme (2c) is an approximation of (lb) . The system (lb) 
is a hyperbolic system of conservation laws with a right hand side which consists 
of a viscosity term in the third equation, a lower order term and a non-local term 
which comes from solving the elliptic equation for $ in terms of n and no [see 9]. 


It is well known that solutions of a system of hyperbolic conservation laws de- 
velop discontinuities after a finite time for general initial data. These discontinuties 
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represent shocks and contact discontinuities. In the present problem, these may 
be slightly smeared because of the (extremely small) viscous terms. There are ad- 
ditional complications due to the stiffness of the source terms. Also these viscous 
shocks can only appear if the Mach number exceeds one. 

There has been a large amount of work in recent years towards developing high 
resolution approximations of the numerical solution to these problems. Schemes of 
the general type (2) are called shock capturing algorithms, since no special treatment 
of the discontinuties is needed. 

First order accurate difference schemes are generally unsuitable for systems 
such as (lb) due to the extreme smearing of the steep gradients in the solution. 
Classical second order accurate schemes, however, are generally unsuitable owing 
to the spurious oscillations which they generate in the neighborhood of discontinu- 
ities. Such oscillations are not only aesthetically undesirable but often lead to a 
breakdown in the stability of the scheme and nonphysical effects such as negative 
densities. 

The feature of the classical schemes which leads to unwanted oscillations is that 
they are based on a constant stencil. Moreover these schemes are linear whenever 
the equations they approximate are linear. In the last fifteen years there has been 
much work on overcoming such problems where the coefficients and the width of 
the stencil depend on the current solution at each time step. Thus the schemes are 
nonlinear even when the equations are linear. This is absolutely necessary [11] for 
the properties of: both high order accuracy in the region of smoothness and the 
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absence of spurious overshoot and undershoots near discontinuities, to be valid. 

The basic method we shall use here which has these desirable properties is 
a modification of the essentially non-oscillatory (ENO) schemes developed in [11]. 
The key ideas are the following: 

In each computational cell the data at t = t n are approximations to the cell 
averages of the vector u, denoted by w". From these cell averages a polynomial 
approximation to u, of order r, is reconstructed in a nonlinear, adaptive fashion. 
This will give us an r-th order accurate method. Values of u™ +3 for s near zero 
are used in such a way as to minimize the likelihood that the stencil crosses over a 
discontinuity. Finally the solution to this piecewise polynomial initial value problem 
is computed approximately, keeping its non-oscillatory properties, and for t n+1 = 
t n + At, At small, the solution is reaveraged over the cell yielding u™ +l . 

As an example, the approximate solution to the scalar equation u t = —u x , 
using a second order accurate ENO scheme is: 

Uj " +1 = «? - A[(«7 + 1*7(1 - A)) - («7_. + 1*7_,(1 - A))] 

where A = ^ < 1 and sj = m[A+u”, A_u"]. for 

X 

m[x,y] = r~\H( x y) min(|x|, |y|) 

l x l 

where H is the Heaviside function . 

The only difference between (2) and a conventional second order accurate ap- 
proximation comes in the definition of the approximate slope s” which is a nonlinear 
function of u" +1 , u”, u”_ x . 
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The scheme we axe vising here is a modified version of the scheme developed in 
[11]. We use the pointwise scheme developed by Shu and Osher [25]. The scheme 
for a scalar one dimensional conservation law is the following. Let the conservation 
law be given by: 

u t + f(u) x = 0 


We partition the x-t space into computational cells of the following type: 




Then any conservative scheme can be written as [15]: 


At 


“f 1 


where the fj_ i represents the numerical flux flowing into the j-th computational 
cell. 

The fluxes can be computed in the following fashion as described in [25]. As- 
suming that there is h(Q such that: 


± f MCK 


then it follows that: 


f(u(x)) x \x=Xj — 


h(x j+ ^)~ h(xj_i 7 ) 


Ax 


To calculate h(x) let: 


H(x) = f h(C)dC 

J — oo 
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Then 


H(x j+i ) = f <+i h(0<K = £ f ' + * MCK = Ax £ /(«(*,)) 

J — OO • ^ */x, 1 




The values of H(x) at the grid points can be calculated from the above and 
then we use a polynomial interpolation with an adaptive stencil to reconstruct the 
H(x) function. We call the reconstructed polynomial Rj + x.(x). Then the flux can 
be computed as: 


•0+* = 

R(x) is constructed in the following way: 


Let us define the following divided differences inductively and using the nota- 
tion that f(j) = f(xj) 

H°(j) = H(x j+i ) 

H\j) = + 1) - H°U)) = f°U) 

H\j) = + 1) - ir-O)) = i/— (j) 

We emphasize that there is no need to calculate the value of H or any of its 
differences since everything can be calculated in terms of the divided differences of 
f(u). To choose the suitable stencil for interpolation we calculate the Roe speed [21] 
according to: 

_ /O' + *) - /O') 

«0 + 1) - «0) 

Let the variable ishift be the amount of shift of the stencil from the j-th grid point; 

if a,j + 1 . > 0 ,then ishift = 0, 


then we compare the divided differences and we move the stencil according to 
if I f\j + ishift) | > | f'(j + ishift — 1)| then ishift=ishift — 1 

We continue in this manner and let m=j-fishift; then the reconstructed poly- 
nomial can be written as: 

R j+ $(*) = H°(m) + H l {m){x - x m _±) + H 2 (m)(x - x m _^)(x - x m+ ^) + • • • 
The value of the flux can be calculated by evaluating the first derivative of R(x) at 

x j+ 1 r 

The above scheme is known to contain nonphysical “expansion shock” solutions 
occasionally. To remedy this case let us define a critical cell where / (u) has a zero 
in the interval / = [uj,u J+1 ]. These are the only cells that could contain expansion 
shocks. 

In a critical cell we use the following scheme: 

Let a j+L = max/ |/ (it)|. Then define the following quantities for j — r + 1 < 
s < j + r where r is the order of the scheme: 

/ + (s) = “(/(«) + a j+ iu(s)) f~(s ) = i(/0) - a >+ ^ti(a)) 

We do the reconstruction for /“ using ishift = 1 and for / + using ishift = 0. 

The flux can be calculated as a sum of the flux moving to the right and the 
flux moving to the left. 

a A 

% ~ fj+ % fj+ % 
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The above scheme can be extended from scalar equations to systems of con- 
servation laws in many ways. We used the following version in our calculations [see 
[25]]. 

Let u t + f(u) x = 0 be the system. Then let represent the Jacobian matrix 
of the system. 

The above system has a complete set of eigenvectors and there is a matrix U 
such that: 

U~ 1 ?fu = A 

ou 

where A is a diagonal matrix. 

Let c 2 = an d 6 = ^ and B = and h = |c 2 + \v 2 . Here, c is the 
velocity of sound waves in the electrons and in this connection we define the Mach 
number of the fluid as M = For our system the above matrices are explicitly: 


111 
U = I m(v — c) mv m(v + c) 
m{—cv + h ) \mv 1 m(cv + h ) 


U~ l = 



— If J_ + if-M 

2 V cm ' m ' 2'm' 

vb 
m 

1/ _1 vb \ 

2 v cm m ) 


m ' 
—6 


1(A) 

2 V m > ■ 


( v — c 0 0 \ 

A = (A ij) = 0 v 0 

\ 0 0 v + c J 


In order to do the field by field decomposition of the fluxes we evaluate the U 
matrix at the Roe average of Uj and u J+1 [see [21]]. 
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Let u.j + 1 = u Roe (uj,Uj+ 1 ) then we calculate 

This gives us the components of the divided differences of the f in the direction 
of the fields;then using the algorithm for scalar equation we calculate the flux for 
each field. We use a J+ i = X(Uj + i) for each field. Then we transform back to the 
original variables via: 

f j+ i = U(u j+k )f j+ 1 

This scheme will be r-th order accurate in space where r is the order of the 
reconstructed polynomial Rj + i.(x). To obtain high order accuracy in time, having 
the ENO property, one may use the Runge-Kutta methods developed in [25]. How- 
ever, since our primary interest here was in steady state calculations we used only 
first order forward differencing in time. The resulting scheme would be Unstable if 
the stencils were fixed', the adaptivity stabilizes the linearly unstable method. 

4. THERMAL EQUILIBRIUM 

If we fix the bias voltage to be zero we expect the system to reach an equilibrium 
state where the temperature of the electron converges to the lattice temperature 
and velocity to zero. 

Under these assumptions the system will be simplified to 

kf,Tn x = en$ r 
= e(n - no) 
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n x (0) = n x (l) = 0 


$( 0 ) = 

e rii 


*(1) = *£ ln( n(!)) 


n i 


Note that we get the following relation: 

n(x) = me** 7 v 

If we let tf(x) = -jAp$(x) and 7 2 = ek fy- then the above system is equivalent to 
the following non-linear boundary value problems: 

f 7 2 ^w - n t -e'* ,(x) = -n D (x) 

l ®*(o) = ®«(i) = o 

f 7 2 ln(n(x)) xx - n(x) = -no(i) 
l n x (0) = n x (l) = 0 

Rather than solving any of the above we choose the following initial conditions and 
let the system evolve in time. 


n(x,0) = n D (x) u(x,0) = 0 T(x,0) = T o 


The calculated solution contains the junction and its electric field. Also, the 
density has to be used as initial condition for the rest of the calculations. The 
numerical results are presented in the numerical section. 

5. PHYSICAL PARAMETERS OF THE DEVICE 

The simulated device is a silicon MOSFET channel. The device is an n + — n— n + 
junction of length one micron. 

We use the following doping profile: 

no( x ) = 5 x 10 17 cm~ 3 if 0 < x < .25 or .75 < x < 1 
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n D (x) = 2 x 10 15 cm 3 


if .35 < x < .65 


The two regions are connected by a properly scaled version of the following poly- 
nomial: 

Q(x ) = — 5x 7 + 21x 5 — 35x 3 + 35x + 16. 


We also use the following conditions at the boundary: 


6 Tt i 

<£(1) = n ^ -) + vbias(t) 

e rii 


n x (0) = n x (l) = 0 Px(0) = p r (l) = 0 


w. 


:(0) = Ml) =0 


We use the following system of units: 

Length : 10 -6 m = 1 micron Time : 10 -12 s = lpico second 

Mass : 10 ~ 30 Kg Potential : Volt Temperature : Kelvin 

Energy : 10 -18 J Charge : 10 -18 C Capacitance : 10 -18 .F 

In the above system of units the constants have the following values: 

m = 0.2 6m e = .26 x .9109 /i 0 = 0.14 e = .1602 v 9 = .1 

k b = .138046 x 10 -4 n, = .014 e = 11.7 x 8.85418 

6. NUMERICAL RESULTS 

In the first experiment we calculate the density of electrons under zero bias. The 
result can be used for calculating the junction capacitance and voltage difference 
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across the junction. Also we use the calculated density as the initial data for the 
rest of the experiments. 

From a numerical point of view this experiment can be used to test the accuracy 
of the scheme. Since the velocity and current must go to zero, the difference can be 
used as a measure of error. 

In our experiment velocity and current converged to a small oscillatory solution 
around the zero value. Also we expect the oscillations in the velocity field to create 
some difficulty for our numerical scheme since it could lead to differencing in the 
wrong direction. 

Indeed, this created spurious oscillations that were eliminated by a small mod- 
ification in the ENO scheme. We restricted the movement of the stencil such that it 
would interpolate across oscillations. This can be done in the following fashion. For 
the sixth order ENO note that the ishift variable takes values between -5 and 1. We 
simply forced it to be between -4 and 0. This procedure would create small oscilla- 
tions across shocks, but one can use a modified version to overcome that difficulty. 
In that version for updating ishift we use: 

if j f*(j 4- ishift) | > |((1 4- e)/*(j 4- ishift — 1)| then ishift =ishift — 1, 

where e is of order of the truncation error of the scheme and positive or negative 
depending on the value of ishift. This would make the choice of stencil biased 
towards the center. In numerical calculations this modification led to a smoother 
solution with fewer spurious oscillations. 
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The results from the the experiment are shown in the fig 2-6. In figure 1 the 
density of donors is shown for reference. As can be seen in figure 2 ,the electrons 
from both sides expand into the channel. This creates electron waves which settle 
down approximately in three pico seconds. The junction is formed and the electric 
field over the junction is created. In figure 4 the steady state of the velocity is 
shown. The amplitude of the oscillation indicates the size of the error in the rest of 
the calculations. 

In the next set of experiments we apply one volt at time zero and calculate the 
evolution of the system in time. The system reaches steady state in approximately 
five pico seconds. In figure 7 the current is shown as a function of time, and in 
figures 8 and '9 the evolution of electric field and velocity are shown. In figures 10 
and 11, for comparison, we duplicate the results obtained by [19], [9]. 

As an approximation to an TV curve we apply a bias voltage in form of a 
“ramp”, Then we graph the calculated current versus voltage in figure 12. If the 
slope of the ramp is small it would be very close to an TV curve. There is some 
oscillation in the curve which is due to the jump in the derivative of the bias voltage. 
If the slope of the ramp is small, it would disappear. In figures 13-17 we show the 
evolution of the solution as the velocity is increased. 

If the temperature of the lattice is raised the sound speed increases and the 
Mach number becomes small. To test the scheme for small Mach number we set 
the lattice temperature at 1000K. The evolution of velocity and temperature are 
shown in the next two graphs. 
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In the next experiment we set the lattice temperature to 50K. This would 

reduce the sound speed . The results are given in figures 20-22. We emphasize that 
* 

we need to modify our choice of stencil to avoid oscillations across steep gradients 
or shocks. 

In the figures 23 and 24 we show just such a situation where the solution 
develops a very sharp profile. The density and the velocity are shown here. The 
shock capturing algorithm is designed to handle such situations. 

In the next experiment we set our bias voltage to be a square wave and then a 
sinusoidal wave. The results are shown in figures 25 -29. 

In the last experiment we compare the calculated solutions in the steady state 
using different orders of accuracy in the x variable. From extensive calculations we 
found out that momentum is the most sensitive variable to error and comparison 
of momentum would show errors that are not obvious in the other variables. 

The calculated results show spurious oscillations for low order schemes that 
(perhaps surprisingly) decrease with increasing order of accuracy. 

7. CONCLUSION 

We presented a numerical scheme that is able to solve the Hydrodynamic model 
efficiently and accurately. 

The scheme was first order in time, but it can be extended to high orders of 
accuracy using the Runge-Kutta methods developed in [25]. Also the high order 
accuracy in space turned out to be essential at the junctions. The equations allow 
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very steep gradients including discontinuous solutions, but there was no observation 
of steady state shocks. However, if a large bias voltage is applied suddenly, shocks 
form which dissappear in a short time. 

There is current work under progress to extend the calculations to holes and 
also to two spatial dimensions. 


FIGURES 


Fig 1: 

density of donors 

Fig 2: 

velocity 


zero bias, Time = 3, ^ = .2, Ax = .01 

Fig 3: 

temperature 


zero bias, Time = 3, ^ = .2, Ax = .01 

Fig 4: 

velocity in steady state 


zero bias, Time = 3, ^ = .2, Ax = .01 

Fig 5: 

n — tid in steady state 


zero bias, T — 5, ££ — .2, Ax = .01 

Fig 6: 

electric field in steady state 


zero bias, Time = 3, ^ = .2, Ax = .01 

Fig 7: 

Current vs. Time 


v(t) = H{t), Time = 5, ££ = .2, Ax = .01 

Fig 8: 

electric field in time 


v(t ) = H(t), Time = 5, = .2, Ax = .01 

Fig 9: 

velocity in time 


22 



v(t) = H(t), Time = 5, ££ = .2, Ax = .01 
Fig 10: comparison of velocity profile for 

bias voltage 1., 1.5, 2. 

Fig 11: comparison of temperature profile for 

bias voltage 1., 1.5, 2. 

Fig 12: Current vs. Voltage 

v(t) = t, Time = 20, ££ = .1, Ax = .01 
Fig 13: Evolution of density 

v(t) = t, Time = 20, = .1, Ax = .01 

Fig 14: Evolution of electric field 

v(t) — t, Time = 20, ^ = .1, Ax = .01 
Fig 15: Evolution of temperature 

v(t) = t, Time = 20, = .1, Ax = .01 

Fig 16: Evolution of velocity 

v(t) = t, Time = 20, ^ = .1, Ax = .01 
Fig 17: Evolution of Mach number 

v(t) = t, Time — 20, = .1, Ax = .01 

Fig 18: Velocity in time 

v(t) = 2 H(t), Time = 3, ££ = .1, Ax = .01, T 0 = 1000A 
Fig 19: evolution of temperature 

v(t) = 2 H(t), Time = 3,j£ = .1, Ax = .01, T 0 = 1000A" 
Fig 20: density in time 
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v(t) = H(t), Time = 4, = .01, Ax = .01, T 0 = 50 K 

Fig 21: evolution of mach number 

v(t ) = H(t), Time = 4, = .01, Ax = .01, T 0 = 50A 

Fig 22: evolution of density 

v(t ) = H(t), Time = 4, ££ = .01, Ax = .01, To = 50 K 
Fig 23: density in time 

v(t) = 2 H(t), Time = 1, = .01, Ax = .01, T 0 = 50 K 

Fig 24: density 

v(t ) = 2 H(t), Time = 1, = .01, Ax = .01, T 0 = 50A 

Fig 25: Current vs. Time 

v(t ) = SQlF(t), Time = 6, ££ = .1, Ax = .01, T 0 = 300A 
Fig 26: velocity in Time 

v(t) = SQlF(t), Time = 6, ££ = .1, Ax = .01, T 0 = 300A 
Fig 27: Current vs. Time 

v(t ) = 2sin(27r<), Time = 6, = .1, Ax = .01, To = 300A” 

Fig 28: density in time 

v(t) = 2sin(27rt), Time = 5, = .1, Ax = .01, To = 300A 

Fig 29: temperature in time 

v(t ) = 2sin(27rt), Time = 5, = .1, Ax = .01, T 0 = 300A 

Fig 30: density, velocity, and momentum not to scale, third order stencil 

Fig 31: density, velocity, and momentum not to scale, sixth ^rder stencil 
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xlO 5 Density of Dnors 



Fig 1: 


density of donors 
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xlO< 



Fig 5: n-no in steady state 

zero bias, T — 5, = .2, Ax = .01 
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volts/micron 


Electric Field 



Fig 6: electric field in steady state 

zero bias, Time = 3, = -2, Ax = .01 
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Current vs. Time 

v(t ) = H(t ), Time = 5, #7 = .2, Ax = 




ELECTRIC FIELD 


Time= 5 Bias Voltage = 1 
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VELOCITY 


Time= 5 Bias Voltage = 1 
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micro amp 



Fig 12: Current vs. Voltage 

v(t) = t, Time = 20, = .1, Ax = .01 
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Fig 15: Evolution of temperature 

v(t) = t. Time = 20, ££ = .1, Ax = .01 
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AJ-IDOTJA 


Time= 3 Bias Voltage = 2 



Fig 18: Velocity in time 

v(t) = 2 H(t), Time = 3, ££ = . 1 , Ax = .01, T 0 = 1000A' 


46 




-DENSITY 


Time= 4 Bias Voltage = 1 



Fig 20: density in time 

»(<) = H{t), Time = 4, = .01, Ax = .01, T 0 = 50K 
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= .01, Ax = .01, To = 50K 
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xlO 5 Evolution of Density 



Fig 22: evolution of density 

v(t) = H(t), Time = 4, ^ = .01, Aar = .01, T 0 = 50A' 
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-DENSITY 


Time= 1 Bias Voltage = 2 



Fig 23: density in time 

v(t) = 2 H(t), Time = 1, g = .01, Ax = .01, T 0 = 50K 
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Fig 24: 


density 

„(t) = 2 tf (t), Time = 1, 


s .01, A j = 01, To — 50^ 
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pico second 


Current vs. Time 

v(t) = SQW(t), Time = 6, = .1, Ax = .01, T 0 = 300# 




VELOCITY 


Time= 6 Bias Voltage = 1 



Fig 26: velocity in Time 

v(t) = SQW(t), Time = 6, g = .1, A* = .01, T 0 = 300# 


54 


n 


,5 


pico second 


Current vs. Time 

u(t) = 2sm(27r<), Time = 6, = .1, Ax = .01, T 0 = 300# 



-DENSITY 


Time= 5 Bias Voltage = 2 



Fig 28: density in time 

v(0 = 2sm(2rrf), Time = 5, ££ = .1, Ax = .01, T 0 = 300 K 
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TEMPERATURE 


Time= 5 Bias Voltage = 2 



Fig 29: temperature in time 


u(0 = 2sin(2irt), Time = 5, = .1, Ax = .01, T 0 = 300K 
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xlO 5 Third order ENO 



Fig 30: density, velocity, and momentum not to scale, third order stencil 
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Sixth order ENO 


xlO 5 



Fig 31: density, velocity, and momentum not to scale, sixth order stencil 
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